function [countobs,countmat]=counter_seldates(funcmod,fullp,y,vecdates,flag_mode)
% 
% Counterfactual for selected dates only 
% Original states fed through all until a certain date, after which only
% either that or all but that shock are on 
%
% Same as KFILTER_ALL but must provide a 2 x 1 vector indicating 
% Beginning of the counterfactual 
% End       of the counterfactual 
% This must be a position in the rows of Y obtained outside from SAMPLE 
% 
% 
% function counter_seldates(funcmod,fullp,Y,flag_allbutone); 

% --------------------------------
% 1. Mode of the counterfactual 
% -------------------------------
% == 1 One shock at a time 
% == 2 All shocks but one 
if flag_mode==1; 
    disp('Counterfactual ONE shock at a time') 
elseif flag_mode==2 
    disp('Counterfactual ALL but ONE shock on') 
else 
    error('FLAG mode must be 1 or 2')
end 

% -------------------------------------
% 2. Check dates of the counterfactual 
% -------------------------------------
nobs=size(y); 
if length(vecdates(:))~=2 
    error('Vecdates must be 2 x 1') 
end 
if vecdates(1) > vecdates(2) 
    error('Vecdates must be increasing') 
end 

[TT,CC,RR,eu,SDX,ZZ]=feval(funcmod,fullp);
if ~isequal(eu,[1;1]);error('Model is not determinate');end
if any(CC~=0)
    disp('Model demeaning the data');
    ydem=y-repmat((ZZ*CC)',[size(y,1) 1]);
else
    ydem=y;
end
[stt,etamat,smooth_st,yfor,vstar,logLnc,countszero]=kfilter_full(ydem,ZZ,TT,RR,SDX);
stt=stt';
[nobs nst]=size(stt);
nx=size(SDX,1); 
disp('Begin Counterfactual')
countmat=zeros(nobs,nst,nx);
countobs=zeros(nobs,size(ZZ,1),nx);
clear temp; 
% Begin Counterfactual
for ii=1:nx
    
    tempeta=etamat';
    if flag_mode==2
        % Shut down that shock only 
        tempeta(ii,vecdates(1):vecdates(2))=zeros(1,vecdates(2)-vecdates(1)+1);
    else
        % Shut down all shocks but that one
        tempind=(1:nx);tempind(ii)=[];
        tempeta(tempind,vecdates(1):vecdates(2))=zeros(nx-1, vecdates(2)-vecdates(1)+1 ) ;
    end
    % ==================================================
    % Loop follows KFILTER_FULL.m
    ysim=zeros(size(ZZ,1),nobs);
    stsim=zeros(nst,nobs);
    stsim(:,1)=smooth_st(1,:);
    ysim(:,1)=ZZ*stsim(:,1);
    for jj=2:nobs;
        stsim(:,jj)=TT*stsim(:,jj-1) + RR*(tempeta(:,jj));
        ysim(:,jj)=ZZ*stsim(:,jj);
    end
    % =====================================================    
    countmat(:,:,ii)=stsim';
    countobs(:,:,ii)=ysim';
end
